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Abstract 



We study hysteresis for a two-dimensional, spin-1/2, nearest-neighbor, ki- 
netic Ising ferromagnet in an oscillating field, using Monte Carlo simulations. 
The period-averaged magnetization is the order parameter for a proposed dy- 
namic phase transition (DPT). To quantify the nature of this transition, we 
present the first finite-size scaling study of the DPT for this model. Evidence 
of a diverging correlation length is given, and we provide estimates of the 
transition frequency and the critical indices f3, 7 and v. 

PACS number(s): 64.60.Ht 64.60.Qb, 75.10.Hk 05.40.+j 



Typeset using REVTgX 



Although nonequilibrium phase transitions have been studied for over two decades, the 
understanding of their universality and scaling properties remains much weaker than for 
equilibrium critical phenomena. An important tool to study equilibrium phase transitions is 
finite-size scaling, which has also been applied to nonequilibrium transitions in which both 
the driving force and the different states are stationary, and the nonequilibrium behavior 
results from the dynamical rules. Examples include diffusive lattice gas models fl]||, Ising 
models with competing kinetic processes at different temperatures P, and systems with 
multiplicative noise |4]]. 

In this letter we present what to our knowledge is the first finite-size scaling analysis of 
dynamic critical phenomena for a system in which the nonequilibrium behavior is due to an 
explicit time-dependence of the Hamiltonian, and for which the nonequilibrium states are 
nonstationary in time. This is the kinetic Ising model in an oscillating field. We focus on 
the nature of the dynamic phase transition (DPT) in this model, which was proposed by 
several workers based on numerical observations @-|7[|. Our results clarify the nature of the 
DPT for this specific model and also provide evidence for the relevance of universality and 
finite-size scaling concepts to dynamic phase transitions in non-stationary systems. 

The relaxation of kinetic Ising models prepared with all spins aligned in a strong field, 
which is suddenly reversed, models the dynamics of metastable phase decay in static field 
HIJ. Such simulations have been used to study magnetization switching in anisotropic 
ferromagnets. The hysteretic response to an oscillating field has been studied by, among 
others, mean-field |10|-p^] and Monte Carlo (MC) [|5| |T|,p~3| ^T^] methods, and some of the 
results have been used to analyze hysteresis loops from experiments on Fe and Co ultrathin 
films fTT|,|18f. The relevance of our results thus extends beyond nonequilibrium statistical 



mechanics to the numerous fields in which kinetic Ising systems are used to model specific 
systems. 

A dynamic phase transition, in which the period-averaged magnetization Q passes from 
a disordered state (Q=0) to an ordered state (Q ^ 0), has been observed in mean- field 
|iOHT2| and MC studies |5|-[7j . The location of this transition depends on temperature, field 
amplitude, and frequency. It can be intuitively understood as a competition between two 
time scales: the period of the external field, 2tc/uj, and the average lifetime, (t(H q )), of the 
metastable phase following instantaneous field reversal from H to — H . If 2tt/uj < (t(H )) 
the magnetization cannot fully switch sign within a single period, and Q 7^ 0. If 2n/uo > 
(t(Ho)) the magnetization follows the field, and Q=0. 

Our previous work on Ising models in sudden field-reversal simulations has identified dis- 
tinct decay regimes in which the decay of the metastable phase proceeds through nucleation 
and growth of one or more compact droplets of the stable phase. These different regimes are 
characterized in great detail in Refs. |||9|]. For this letter, the most important result from 
past work is that the metastable decay mode may change drastically as the temperature, 
field strength, and/or system size are varied. At weak fields and/or small system size, the 
decay proceeds by the nucleation and growth of a single droplet of the stable phase [single- 
droplet (SD) regime]. For stronger fields and/or larger systems many droplets contribute 
to the metastable decay [mult i- droplet (MD) regime]. The MD decay mode can be accu- 
rately described by the classical 'Avrami's law" for nucleation and growth [JL9]. For each of 
these regimes, the statistical properties of the lifetime of the metastable phase are different, 
leading to very different responses in both static and time-varying [p||n| fields. 
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Here we emphasize two main points. First, care must be exercised in determining the 
temperature and field dependence of the DPT boundary for a fixed frequency, because the 
response of the model can change qualitatively when the temperature and field strength are 
changed. In fact, our simulations show evidence of a DPT for the MD regime only. Second, 
if the DPT in the kinetic Ising model is a true critical phenomenon, then all of the theoretical 
machinery of finite- size scaling for equilibrium phase transitions should be applicable. 

The model used is a kinetic, nearest-neighbor Ising ferromagnet on a square lattice with 
periodic boundary conditions and Hamiltonian Ti.=—JJ2(ij) s i s j ~ H(t) Ya Si- Here Sj=±l, 
Y{ij) runs over all nearest-neighbor pairs, and J2i runs over all L 2 lattice sites. The ferro- 
magnetic coupling is J—l. Each spin is subject to an oscillating field H(t)=—H sm(ujt). 
We measure the time-dependent magnetization per site, m(t) = (l/ L 2 )Y,i=i Si(t), using 
the Glauber single- spin-flip Monte Carlo dynamic |2(J with updates at randomly cho- 



sen sites. Each attempted spin flip from s, to — s, is accepted with probability W(si — > 
— Sj)=exp(— /3AEi)/[(l + exp(— /3A£?j))]. Here AEi is the energy change of the system re- 
sulting from an accepted spin flip, and (3 = where fee is Boltzmann's constant. The 
time unit is one Monte Carlo step per spin (MCSS). 

All simulations are performed for three system sizes: L=64, 90, and 128 at T=0.8T C . 
This temperature is sufficiently far below T c , the critical temperature of the Ising model 
in zero field, that the thermal correlation length is small compared to L and the critical 
droplet radius. We choose H =0.3J. This field amplitude is such that, for simulations at 
0.8T C in a static field of H=Hq, all three system sizes are in the MD regime. To obtain the 
raw time-series data, the system was initially prepared with either a random arrangement 
of up and down spins with m(t = 0) ~ 0, or with a uniform arrangement with all spins up, 
m(t = 0)=1. We recorded m(t) for several values of u for approximately 1.7 x 10 7 MCSS 
(for the lowest frequencies, 5.9 x 10 5 MCSS). Each of these raw data files store t, H(t), and 
m(t) in increments of 1 MCSS and contain thousands of field periods. The largest files are 
about 800 megabytes and required 9 days (one month) to run for L=64 (L=128) on a single 
node of an IBM sp2. This is one of the most extensive MC simulations of hysteresis in Ising 
systems to date. 

It is useful to think of hysteresis as a competition between two time scales: 2tt/uj and 
(t(Hq)). Therefore we report all our results in terms of the dimensionless period, 

R = (2tt/u;)/ (r(ifo)) • (1) 

The average lifetime has been measured in field reversal simulations to be (t(H )) « 74.5 
MCSS and is independent of L in the MD regime. 

The order parameter of the dynamic phase transition is the period-averaged magnetiza- 
tion, 

M —/';;;,/,,//. (2) 

This definition removes the field oscillations, so that Q is a stationary process. For each 
frequency we obtain the probability density of Q by constructing a histogram of the Q values 
calculated from each period in the corresponding time series. 

Figure [l] shows the probability densities of Q, P(Q), for all three system sizes at a 
frequency near the transition. (The details of locating the transition frequency are given 
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below.) Correlation times that are significant portions of the total run length are manifest 
in the remaining asymmetry of the distributions. Estimating the correlation time from the 
asymmetry in P(Q) we find it to be between 1 and 10 percent of the total simulation length 
near the transition. Away from the transition it decreases rapidly. The behavior is remi- 
niscent of the critical slowing-down seen in equilibrium simulations, and even our extensive 
simulations are not long enough to be fully ergodic. The asymmetry in the distributions is 
not sensitive to the initial condition of the time series. 

At a second-order phase transition there is a divergence in the susceptibility. For equilib- 
rium systems, the fluctuation-dissipation theorem relates the susceptibility to fluctuations 
in the order parameter. For the present system, it is not obvious what the field conjugate 
to Q might be. Therefore we cannot measure the susceptibility directly. However, we can 
calculate the variance in \Q\ as a function of frequency and study its system size dependence. 
We define X as 



X = L Var(|Q|) = L [(Q ) - (\Q\Y\ • (3) 

If the system obeyed a fluctuation-dissipation relation, X would be proportional to the 
susceptibility and both would scale with L in the same manner. Figure |2] shows X vs 1/R 
for three system sizes. For all three values of L, X displays a prominent peak which increases 
in height with increasing system size. This clearly shows finite-size effects in X and implies 
the existence of a divergent length associated with the order-parameter correlation function 
near the dynamic transition. The observation that P(|Q|) displays no peak near |Q|=0 in 
the ordered phase is additional evidence for the second-order (as opposed to first-order) 



nature of this transition [pi 



One would like to find the critical exponents associated with this transition, as well as 
the frequency at which the transition occurs. For the Ising model in zero field, T c is exactly 
known. Then, one can use scaling relations which depend on T c to directly calculate the 
critical exponents. In the present case, no exact solution exists for the transition frequency, 
and the scaling relations for X involve the peak heights and positions. Both of these are 
difficult to measure accurately, even from our extensive data. The cumulant intersection 



method [20,22] is useful for determining the location of a second-order transition when the 



critical exponents are not known. We define the "dynamic" fourth-order cumulant ratio as 

U L = 1-M\. (4) 

H\Q\ 2 )l 

where (\Q\ n )=f™ \Q\ n P{\Q\)d\Q\. Figure | shows U L vs 1/R. Above the transition fre- 
quency, in the (\Q\) ^ phase, Ul approaches 2/3. Below the transition frequency, in the 
(Q)=0 phase, Ul approaches 0. At the transition, the cumulant should have a non-trivial, 
fixed value, U*. Therefore, the location of the cumulant intersection gives an estimate of 
the transition frequency without foreknowledge of the critical exponents. Due to the large 
spacing of our data and possible correction-to-scaling effects, we cannot identify a unique 
intersection point. We estimate the location of the intersection by the crossing of the two 
largest system sizes near 1/R C ~ 0.2910 with U L =U* « 0.61. 

Having estimated the transition frequency, -R" 1 , we can approximate the critical expo- 
nents [ ^0| , p^1 . We obtain estimates for (3/u, 7/V, and \/v using the two largest system 
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sizes. At the transition we use scaling relations for the moments of the order parameter 
((|<5| n ) oc L~ n ^l v ^\ the maximum value of the order parameter fluctuations (A max oc L 7 ^ v ), 
and the position Rj} of the maximum value of the order parameter fluctuations for a partic- 
ular system size (l-R^ 1 — R^\ oc L~ l l u ). Using the scaling relations for either the second or 
fourth moments of \Q\ we estimate (/3/V) n=2 ~ 0.111 and (/3/z/) n=4 « 0.113. Our estimates 
for the other exponents are y/u ~ 1.84 and u ~ 1.1. Simulations with larger system sizes 
would be computationally prohibitive, and smaller system sizes would no longer be in the 
MD regime. 

The scaled probability distributions of \Q\ are given in the inset of Fig. [I] after symmetriz- 
ing and scaling to demonstrate data collapse. The symmetrizing is equivalent to calculating 
the distribution for \Q\. The scaling form is derived in a fashion analogous to equilibrium 



finite-size scaling analysis of order-parameter distributions p0| , |22| . At the transition, we 
assume that the mean of the order parameter scales with L and define the scaling variable 
Q=Ll 3 l v \Q\. Hence, the scaled probability density for \Q\ is given by 

P L {Q) = L-^P{\Q\) . (5) 

The peak positions scale fairly well, the peak heights less so. This could be due to the 
following reasons. The frequency might be sufficiently far from the transition that single- 
parameter scaling is not adequate. There might be corrections to the finite-size scaling that 
are large for these relatively small system sizes. Also, the lack of scaling for the peak heights 
could be due to the asymmetry in P(Q) near the transition. Much longer simulations on 
larger lattices would be needed to resolve this issue. 

We have also carried out an extensive study of hysteresis in a kinetic Ising model in the 
SD regime where we have found evidence of stochastic resonance, but no sign of a dynamic 



phase transition [[L3[]. In the introduction we emphasized that the crossover between the SD 
and MD decay regimes depends on temperature, field strength, and system size. Therefore, 
the very existence of the dynamic transition, as well as the details of its critical behavior, 
may depend sensitively on all of these parameters. 

In conclusion, we have performed a finite-size scaling study of a kinetic Ising model in a 
sinusoidal field in order to clarify the nature of the dynamic phase transition conjectured by 
several authors |^[l3]|nj . For this letter, we emphasize that all simulations were performed 



in the MD regime. The behavior of the order-parameter fluctuation, X, suggests a divergent 
correlation length near the transition frequency. This behavior motivates the application 
of finite-size scaling techniques for second-order phase transitions, analogous to those used 
to describe the ferromagnetic/paramagnetic transition in the Ising model in zero field. We 
use the cumulant of the order-parameter distributions to estimate the value of the transi- 
tion frequency, 1/R C ~ 0.2910. Using scaling relations for the moments and fluctuations 
of the order parameter we estimate the critical exponents to be (3/u ~ 0.11, y/u « 1.84 
and v ~ 1.1. Our results for (3/u and y/u are close to the two-dimensional Ising values for 
the analogous exponent ratios. The result, 2{(3/u) + (y/u) ~ 2.06 ~ d, gives a tantalizing 
indication that hyperscaling may be obeyed. Also, our value for the cumulant intersection, 
U* ~ 0.61, agrees with an extremely precise transfer matrix calculation of U*=0. 6106901(5) 
for the two-dimensional Ising model ||23|| . To precisely calculate the critical frequency, critical 
exponents, and determine the universality class of the transition would require simulations 
for more frequencies, larger systems, and for immensely longer run times to improve the 
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statistics. More detailed analysis of this problem could resolve specific questions about the 
DPT for the Ising model, such as the existence of a tricritical point in the dynamic phase di- 
agram, in addition to elucidating general questions concerning the nature of nonequilibrium 
phase transitions in models with explicitly time-dependent Hamiltonians. 

We would like to thank G. Brown, W. Janke, W. Klein, M. Kolesik, G. Korniss, M. 
Acharyya, and J. Vihals for useful discussions. Supported in part by the FSU Center for 
Materials Research and Technology (MARTECH), by the FSU Supercomputer Computa- 
tions Research Institute (SCRI) under DOE Contract No. DE-FC05-85ER25000, and by 
NSF Grants No. DMR-9315969, DMR-9634873, and DMR-9520325. 
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FIGURES 
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FIG. 1. Probability densities of the period-averaged magnetization, Q=(u /2ir) § m(t)dt for 
three system sizes, L=64 (squares), 90 (triangles), and 128 (stars). The frequency of the field, 
l/i?=0.2910 is near the transition. The asymmetric distributions indicate correlation times on the 
order of the simulation time even for our extremely long runs of 1.7 x 10 MCSS. Inset: Scaling 
function L-M V P(\Q\) vs lPl v \Q\. The value of the scaling exponent is {j3/u) n=2 ~ 0.111. 
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FIG. 2. L 2 Var|Q| vs dimensionless frequency, 1/R. The "disordered phase," ((|Q|)=0), lies on 
the low-frequency side of the peaks. The "ordered phase", ((|Q|) 7^ 0), lies on the high-frequency 
side. Lines connecting data points are guides to the eye. The statistical error bars are estimated 
by partitioning the data into ten blocks. Error bars smaller than the symbol sizes are not shown. 
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FIG. 3. Fourth-order cumulant ratio Ul vs dimensionless frequency, 1/R. for L=64, 90, and 
128. We use the same system size symbols as in Fig. |]. The horizontal line marks Ul=2/3. Lines 
connecting the data points are guides to the eye. Inset: area close to the cumulant crossing. 
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